BZOJ4126 国王奇遇记加强版之再加强版 <多项式插值>

Problem

国王奇遇记加强版之再加强版

Description

Input

共一行，包括两个正整数和。

Output

共一行，为所求表达式的值对取模的值。

Sample Input

5 3

Sample Output

Hint

，

Source

标签：多项式插值

Solution

好题，看了中的题解才懂。以下题解全部部分摘自特殊多项式在整点上的线性插值方法和BZOJ-3157. 国王奇遇记。

1. 多项式整点插值

观察二项式系数，其为一个的次多项式。对于，由于其次数互不相同，故其线性无关。可以发现这个多项式是次多项式线性空间的一组基。
于是对于，其一定可以表示为这个多项式的线性组合，即

由于时，，于是当时，有

根据二项式定理，可知，应用其对进行二项式反演，得

将带入，得

化简一下二项式系数

于是

即

将后面的部分进一步化简，即化简形如的式子。
根据上指标反转，可将其化简

带入得

将后面的再次上指标反转得

即

故对于次多项式，当时，如果能得到的值，可以求出的值。
具体地，两个二项式系数的乘积为

对于分式部分，可以预处理阶乘及其逆元以计算，对于两个乘式，可以预处理的前缀积和后缀积以计算。总时间复杂度为。

2. 幂和

记，求出较小时的通项，瞎猜发现当时，一定有如下形式：

其中是一个次多项式。根据前面多项式整点线性插值，只需求出即可求得。
注意到在中被消去，于是可以得到的递推式：

于是可以将表示为的形式，然而还无法求出。
根据多项式插值时推导的公式，当时，有

将作为带入得

将移到右边得

将带入即可解出，复杂度。然后通过前面插值可求出，带入即可求得。总时间复杂度。
注意由于时不满足性质，需要单独计算；此外还需暴力计算的情况。

简单版之再简单版见BZOJ3157，简单版见BZOJ3516。

Code

#include <bits/stdc++.h>
#define MAX_M 500000
#define P 1000000007
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int n, m; lnt F[MAX_M+5];
lnt pw[MAX_M+5], fac[MAX_M+5], inv[MAX_M+5];
bool NotPri[MAX_M+5]; vector <int> pri;
lnt Pow(lnt x, lnt k) {
    lnt ret = 1;
    for (; k; k >>= 1, x = x*x%P)
        if (k&1) ret = ret*x%P;
    return ret;
}
void init(int n) {
    fac[0] = inv[0] = inv[1] = pw[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!NotPri[i]) pri.push_back(i), pw[i] = Pow(i, m);
        for (int j = 0; j < (int)pri.size(); j++) {
            if (i*pri[j] > n) break;
            NotPri[i*pri[j]] = true;
            pw[i*pri[j]] = pw[i]*pw[pri[j]]%P;
            if (!(i%pri[j])) break;
        }
    }
    for (int i = 1; i <= n; i++) fac[i] = fac[i-1]*i%P;
    for (int i = 2; i <= n; i++) inv[i] = (P-P/i*inv[P%i]%P)%P;
    for (int i = 2; i <= n; i++) inv[i] = inv[i-1]*inv[i]%P;
}
lnt C(int n, int m) {return fac[n]*inv[m]%P*inv[n-m]%P;}
void get_Pnt_Val() {
    lnt k[MAX_M+5], b[MAX_M+5], invm;
    k[0] = 1, b[0] = 0, invm = Pow(m, P-2);
    for (int i = 1; i <= m+1; i++)
        k[i] = k[i-1]*invm%P, 
        b[i] = (b[i-1]*invm%P+pw[i])%P;
    int f = (m&1) ? -1 : 1; k[0] *= f, f = -f;
    for (int i = 1; i <= m+1; i++, f = -f)
        k[0] = (k[0]+f*C(m+1, i)*k[i]%P)%P, 
        b[0] = (b[0]+f*C(m+1, i)*b[i]%P)%P;
    F[0] = -b[0]*Pow(k[0], P-2)%P;
    for (int i = 1; i <= m; i++)
        F[i] = (k[i]*F[0]%P+b[i])%P;
}
#define getL(i) (i < m ? pre[m-i-1] : 1)
#define getR(i) (i > 0 ? suc[m-i+1] : 1)
lnt Poly_Inter() {
    lnt ret = 0;
    lnt pre[MAX_M+5], suc[MAX_M+5]; pre[0] = n-m, suc[m] = n;
    for (int i = 1; i <= m; i++) pre[i] = pre[i-1]*(n-m+i)%P;
    for (int i = m-1; i >= 1; i--) suc[i] = suc[i+1]*(n-m+i)%P;
    for (int i = 0, f = (m&1) ? -1 : 1; i <= m; i++, f = -f)
        ret = (ret+f*getL(i)*getR(i)%P*inv[i]%P*inv[m-i]%P*F[i]%P)%P;
    return ((Pow(m, n)*ret%P-F[0])%P+P)%P;
}
int main() {
    read(n), read(m), init(m+1);
    if (m == 1)
        printf("%lld\n", 1LL*n*(n+1)/2%P);
    else if (n <= m) {
        lnt pwm = m, sum = 0;
        for (int i = 1; i <= n; i++, pwm = pwm*m%P)
            sum = (sum+pw[i]*pwm%P)%P;
        printf("%lld\n", sum);
    } else
        get_Pnt_Val(), printf("%lld\n", Poly_Inter());
    return 0;
}